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G\ ' ABSTRACT 

, We present interferometric observations of N2H+(l-0) in the starless, dense 

I core L1544 in Taurus. Red-shifted self-absorption, indicative of inward motions. 
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is found toward the center of an elongated core. The data are fit by a non- 
spherical model consisting of two isothermal, rotating, centrally condensed layers. 
Through a hybrid global-individual fit to the spectra, we map the variation of 



Tij- I infall speed at scales ~ 1400 AU and find values ~ 0.08 km s ^ around the core 

' center. The inward motions are small in comparison to thermal, rotational, and 



gravitational speeds but are large enough to suggest that L1544 is very close to 



g ■ forming a sta.. 
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I ' Subject headings: ISM: individual(L1544) — ISM: kinematics and dynamics 

i-t • stars: formation 



1. Introduction 

The very first stage in models of the star formation process is the gravitational collapse of 
a starless dense molecular core (e.g., Shu, Adams, & Lizano 1987) but this stage is relatively 
short-lived, and therefore rare. One of the few known examples is the dense core, L1544, in 
Taurus which shows no evidence of an embedded young stellar object but appears to have 
extended inward motions (Tafalla et al. 1998; hereafter T98). This core is therefore an 
excellent testing ground for theories of low mass, isolated, star formation. 

T98 presented single-dish molecular line observations and analyzed the large scale 
structure and dynamics of L1544. They find that lines of CO and its isotopes are single 
peaked and trace an outer envelope, but higher- dip ole- moment lines of H2CO, CS, and N2H"'" 
are strongly self-absorbed and trace a central core. The self-absorption is predominantly red- 
shifted as expected for inward motions (e.g., Myers et al. 1996; hereafter M96) although there 
is also a small region of blue-shifted reversals toward the south. The inferred inward speeds 
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range from 0.10 km in the northeast to —0.03 km (relative outward motions) in the 
region of reversals in the south. T98 note that the dynamics differ from several theoretical 
predictions for core collapse: the infall speeds are too large for ambipolar diffusion (Ciolek 
& Mouschovias 1995); the spatial extent of infall is too large for inside-out collapse (Shu 
1977) at such an early stage; and the density is too highly centrally condensed to be due to 
collapse from a uniform density configuration (Larson 1969). 

In order to measure the structure and dynamics in greater detail, we made high 
resolution observations of N2H"'"(l-0). This line was chosen because the T98 observations 
show a compact, centrally condensed core with double-peaked spectra. These observations 
are the first kinematic measurements of a starless core on the size scale ~ 2000 AU 
(~ 0.01 pc). Our main finding is that the maximum infall speed is (surprisingly) essentially 
the same, ~ 0.1 km s~^, as on much larger scales. 

2. Observations 

L1544 was observed with the ten-antenna Berkeley-Illinois-Maryland array [| (BIMA) in 
its compact C configuration for two 8 hour tracks on November 4 and 5, 1997. Projected 
baselines ranged from 1.9 to 25.1 kA. The phase center was a(2000) = 5^04"^16!62, 5(2000) = 
25°10'47'.'8. Amplitude and phase were calibrated using 3 minute observations of 0530-1-135 
interleaved with each 25 minute integration on source. The passband shape was determined 
from observing the bright quasar 3C 454.3. No suitable planets were available for flux 
calibration and we therefore assumed a flux density of 2.4 Jy for 0530-1-135 based on flux 
density measurements taken within a month of these observations. From the scatter in these 
measurements, the flux density scale is accurate to 30%. 

The digital correlator was configured with 512 channel windows at a bandwidth 
of 6.25 MHz (0.04 km s~^ velocity resolution per channel) centered on the 7 hyperfine 
components of N2H"''(l-0) (Caselli, Myers, & Thaddeus 1995) in the lower sideband and 
C^^S(2-1) in the upper sideband. Eight 32 channel windows at a bandwidth of 100 MHz. 
were used to measure continuum radiation. Data were reduced with the MIRIAD package 
using standard procedures. The data sets from each day were calibrated separately and 
then transformed together using natural weighting. The resulting "dirty" map was cleaned 
and restored with a gaussian beam of FWHM size 14'.'8 x 6'.'6 at position angle —1°. The 
continuum map shows a 5 mJy peak, corresponding to 5 o", at offset position, Aa = 4", A5 = 
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12", but this has not been verified by the IRAM Plateau de Bure interferometer (J.-F. Panis, 
personal communication) and may be due to contamination of low level line emission from 
species such as CH3OH in the wide windows. We also searched for continuum emission at 
3.6 cm from a very young protostellar wind using the VLA in B-configuration on July 5, 
1998. No emission was detected to a 3o" level of 39 fiJj which is considerably less than the 
outflows detected by Anglada (1995) and, extrapolating his data, implies an upper limit to 
the bolometric luminosity of any protostar in L1544 of 0.03 Lq. 

The N2H+ lines were detected with high signal-to-noise but no emission was detected in 
the other spectral window implying a 3 cr upper limit of 0.25 Kkms^^ for the C^^S emission 
integrated over velocities 6.8 to 7.6 km s~^. A map of the velocity integrated N2H+ emission 
for the isolated hyperfine component FiF= 01 — 12 is displayed in Figure 1. There is a 
single, compact, elongated structure with FWHM size 7000 AUx3000 AU. Double-peaked 
spectra are found toward the map center as in T98, but before analyzing these profiles it is 
necessary to account for the effect of the missing zero-spacing flux (see Gueth et al. 1997). 
Therefore, we combined the single-dish IRAM 30 m data of T98 with the BIMA map: the 
single-dish data were shifted by 13" to the south and then scaled to match the interferometer 
data in the region of visibility overlap (6 m to 30 m) using a best fit gain 4.8 Jy K^^. The 
resolution of the resulting data set is the same as for the interferometer data alone but a 
circular restoring beam of FWHM size 10" (having approximately the same beam area) was 
used to create a regular grid of spectra for analysis. 



3. Analysis 

The combined IRAM-BIMA spectra are also double-peaked toward the map center and 
we attribute this to self- absorption for the reasons outlined in T98. The red peak is stronger 
than the blue as expected for a contracting core but the contrast between the two is not as 
great as in the interferometer map alone which implies that the absorbing region is more 
spatially extended than the emitting region. 

To examine the structure and dynamics of the core in further detail, we fit the data 
with a simple model of a collapsing core. The self- absorption implies a decreasing excitation 
temperature gradient away from the core center toward the observer. Many previous 
models of core collapse (e.g., Zhou 1995) have assumed a spherical gas distribution, but 
Figure 1 shows that the core is highly elongated with an axial ratio ~ 0.4, so such models 
are inappropriate in this case. Here, we allow for variations in parameters parallel and 
perpendicular to a major axis in the plane of the sky, and consider two isothermal layers, each 
of constant thickness, but at different densities so as to provide the excitation temperature 
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inversion. We were motivated to consider this heuristic model because it greatly simplifies 
the treatment of the radiative transfer yet fits individual spectra quite well. It allows us 
to infer the variation of density in the plane of the sky and the magnitude and spatial 
distribution of the relative motion between the two layers (the infall speed) . 

The fine brightness temperature is given by M96 equation (2) and depends on the 
excitation temperature and optical depth of each layer. The excitation temperatures are 
derived using a two level approximation as in M96 equation (8a,b) and the peak optical 
depth of the 101-012 transition is 



for Ti = rpk,iexp[— — Vi)^/2(j^], where Vi is the systemic velocity and a is the velocity 
dispersion. Here, the index z = / or r for the front or rear layer respectively, and X is the 
abundance of N2H+ relative to H2, rij is the density of H2, A.Zi the layer thickness, and Tex,i 
the excitation temperature. 

This formalism allows each spectrum to be fit very well but requires at least 6 
independent parameters for each spectrum (M96). Instead, we present a global fit to the 
data that describes the core structure and dynamics with the addition of fewer parameters. 
The route to a global fit was guided by the individual fits; these showed that the velocity 
dispersion did not vary greatly but the excitation temperature and optical depth both 
increased sharply toward the map peak. The kinetic temperature is constrained to the 
range Tk = 11 — 13.5 K by the CO observations of T98 and is insufficient to account for the 
range of excitation temperature, implying that the density must increase toward the map 
center. We therefore fixed Tk and a and chose the functional form. 



to describe the spatial variation of the density, ni{x,y) = ripi^^i/ f{x,y), where x — y defines 
an orthogonal coordinate system rotated with respect to the a — S observational frame by 
an angle 9 (measured anti-clockwise from north), chosen so that Ax > Ay (i.e., x is the 
major axis). The kinetic temperature and velocity dispersion are then held constant, both 
spatially within a layer and from layer to layer. 

The individual fits (and also channel maps) show that there is a significant velocity 
gradient that is approximately uniform and along the major axis and therefore we set 
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Vr{x) =Vo + {dv/dx){x - Xpk), Vf{x,y) = Vr{x) + Vir,{x,y), 
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where Vq is a constant velocity offset, dv/dx is a constant velocity gradient {dv/dy = 0), and 
Vin — i^f—Vr is the velocity difference between the front and rear layers (the infall speed). 

In addition to Vi^, there are, formally, a total of 14 parameters in the final model: 
x-pk-iU-pk-iO to define the coordinate system; npk,f, Az/, npk,r, Az^, Ax, Ay to describe the 
densities in each layer; the velocity parameters, Vq and dv/dx; and the constants, T]^,a, and 
X. This is far fewer than the >200 parameters needed to fit the spectra in Figure 2 using 6 
parameters per spectrum. Furthermore, the parameter space is quite tightly constrained by 
the maps showing the approximate center, angle, and velocity gradient of the core, and the 
estimates of kinetic temperature, sizes, and linewidths that were determined by T98. Also, 
not all the parameters are independent; all values of Azi and X with the same product result 
in the same model output. The principal parameters that were varied were the six density 
parameters and the infall speed. 

We proceeded by assuming a constant value for Vin{x, y), and searched for the best global 
fit by least squares minimization of the difference between the model and observed spectra 
for a grid of 5 x 7 spectra at 10" spacing about the core center (Figure 2). This showed that 
any viable fit must satisfy the following conditions: (1) npk,r ^ Ud ^ npk,f and Az^ ^ Azf 
where ricr is the critical density; the excitation temperature must be high in the rear which 
emits strongly, and low in the front which emits weakly but has a large line-of-sight thickness 
and absorbs strongly. (2) Ax, Ay < (core size); the excitation temperature and optical depth 
must decrease rapidly away from the map center. (3) dv / dx ^ 0; a constant systemic velocity 
is a very poor fit and a velocity gradient is required. (4) fin > 0; the front layer must be 
moving toward the rear layer to account for the asymmetry of the self-absorption. In fact, 
many of the spectra do not show two distinct peaks but rather a bright red peak and a 
blue shoulder indicating that the velocity of the absorption must be shifted relative to the 
emission by an amount approximately equal to the velocity dispersion, i>in ~ a. 

Model spectra are plotted with the observations in Figure 2 for the best fit parameter 
values listed in Table 1. We chose a high value for the abundance, X — Ix 10~^, so that the 
front layer (with peak density npk,f = 10^ cm~^) has a symmetrical size, 2A2;/ = 0.4 pc, that 
approximately matches the region of C^^O emission in T98. The infall speed was determined 
for each spectrum individually by least squares minimization after the other parameters were 
fixed. This hybrid global-individual approach was chosen to emphasize the spatial variation 
of Win. However, its determination requires high signal-to-noise ratios and optical depth, and 
it could only be reliably measured around the map center where there are bright lines and 
strong self-absorption. We find that the best fit has fin increasing toward the core center 
and along the major axis, but because the optical depth falls off rapidly, the significance of 
this fit is only marginally better than a fit with constant Vin — 0.075 km s~^. The difference 
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between the model spectra and the data has a mean of zero and standard deviation 0.12 K 
per 0.04 km channel, about 50% greater than the rms noise in the spectra. The spatial 
variation of column density and infall speed are displayed along with the spectra in Figure 2. 

4. Discussion 

The observations reported here are the first to define the inward motions in a starless 
core on the size scale of star formation, i.e. on a size scale which encloses about a stellar 
mass. We therefore comment here on the model parameters, the inferred core motions and 
their physical basis, and the evolutionary status of L1544. 

Our model properties agree with four independent constraints. We find agreement 
between the size scale, Azj, as discussed above, and the kinetic temperature, Tk = 12 K, from 
the T98 CO observations. The core mass can be derived by integrating the density profile 
out to a distance equal to the spatial thickness of each layer: Mf = 0.2 Mq, Mr = 0.4 Mq 
implying a total symmetric mass 2{Mf + Mr) = 1.2 Mq, which is comparable to estimates 
of the dense gas mass discussed in T98. The density profile and peak column density are 
consistent with dust continuum measurements of the structure of a similar prestellar core, 
L1689B, in the p Ophiuchus cloud (Andre, Ward-Thompson, & Motte 1996). We have also 
fit the data with the exponent in equation (2) varying from 1.8 to 2.5 (with corresponding 
changes in Ax and Ay) without significantly increasing the least squares difference. 

The core dynamics are of greatest interest. The dispersion, a = (a^ + cr|frp)^/^, is 
the quadrature sum of a thermal component, — (/cTk//^)^/^, where k is the Boltzmann 
constant and /i is the molecular mass, and a non-thermal component, unt, which has a 
best fit value of 0.085 km s-^ similar to that measured for C^^^'S by T98. For Tk = 12 K, 
aT(B.2) = 0.22 km s~^, much greater than ctnt and comparable to gravitational speeds 
~ {GMr/ AzrY^"^ = 0.27 km s~^. Therefore, the turbulent pressure support is negligible. 
However, the non-spherical core shape implies that it cannot be supported entirely by 
thermal pressure and there must also be an anisotropic force resisting self-gravity. The 
velocity gradient along the major axis, dv/dx = 3.8 km s~^ pc^\ implies a change in velocity 
Avrot = 0.13/sim km s"^ over the 7000 AU major axis FWHM core diameter, where i is 
the inclination of the gradient to our line of sight. If the non-sphericity is due to rotation, 
however, the aspect ratio constrains cosi < 0.4, so Avrot < 0.14 km s~^ <^ ct(H2) and 
we conclude that rotation is dynamically insignificant. This leaves magnetic fields as the 
most likely explanation for the extended core shape. The observed aspect ratio implies 
approximate equipartition between thermal and static magnetic pressure support (Li & Shu 
1996), corresponding to a field strength 5 ~ 30 /iG at number densities ~ 10^ cm~^. 
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The thermal, gravitational, and rotational speeds are all greater than the maximum 
infall speed. We have determined a range, ~ 0.02 to 0.09 km s""*^, that is similar to the 
values derived from the CS/C^^S analysis by T98 even though the maps are quite different in 
spatial scale. The interferometer map presented here has finer resolution, but covers a much 
smaller spatial extent than the T98 CS map. Comparing the two maps, the high resolution 
observations indicate a sharp increase in density, consistent with free-fall collapse, but no 
corresponding increase in infall speed. 

A possible explanation is that starless cores, such as L1544, are predicted to evolve 
through a scries of quasi-static equilibria (see discussion in Li & Shu 1996). For example, 
neutrals in an outer envelope diffuse through magnetic field lines onto a central core which 
is supported by thermal pressure. The piling up of material on the central core causes its 
density to increase and its (Jeans) size to decrease. The mass-to-flux ratio increases until 
the collapse becomes supercritical and magnetic field lines are dragged into the core which 
remains (quasi-)static. The infall velocities remain sub-sonic except for a small region at 
the envelope-static core boundary at very late times (Li 1998; Basu & Mouschovias 1994). 
As with T98, however, the quantitative comparison shows discrepancies. For example, both 
Ciolek & Mouschovias (1995) and Basu & Mouschovias (1994) predict smaller infall velocities 
than measured here, i.e., ^ 0.02 km s~^, at size scales of 0.02 pc and densities of 4x 10^ cm~^. 
The large scale (~ 0.1 pc) inward motions observed by T98, possibly driven by turbulent 
pressure gradients (Myers & Lazarian 1998) and not included in the above models, probably 
increase the infall speed at all size and density scales. 

N2H"'" is an ion and therefore moves only along, or with, the magnetic field lines. Since 
the ions and neutrals are moving together at a speed much greater than their relative 
drift speed, either we are looking down along the field lines, or the magnetic field lines 
are being dragged in along with the gas. The geometry suggests that the field lines he more 
closely perpendicular than parallel to our hne of sight, and therefore the second possibility, 
supercritical collapse, seems more reasonable. In this case, L1544 may be very close to 
forming a star: Li (1998) shows that, for a spherical cloud, both the ion and neutral infall 
speeds are very small for times t <^ 15rg, where rg ^ 0.4 Myr is the free-fall timescale, but 
then rapidly increase for t ^ Idrg and the core forms a star by 20rflf. There may be many 
cores that have begun collapsing in Taurus, but these calculations suggest that only those 
that are more than 95% of the way there, such as L1544 appears to be, may show detectable 
inward motions. Observational verification, however, awaits the statistics from surveys for 
inward motions in starless cores. 

Discussions with Zhi-Yun Li and Shantanu Basu are gratefully acknowledged. This 
research was partially supported by NASA Origins grant NAGW-3401. 
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Fig. 1. — BIMA map of the integrated intensity map of N2H"'"(101-012) emission in L1544. 
The emission has been summed over velocity channels from 6.8 to 7.6 km s^^. Contours 
begin at, and are in increments of, 0.12 K km s^^ (the dotted contour is at -0.12 K km s^^). 
The beam size is 14'.'8 x 6'.'6 at a position angle of —1°, and is indicated in the lower right 
corner. 




Fig. 2. — Spectra and model fits for tlie combined BIMA-IRAM map. Coordinates are 
shown offset from a(2000) = 5*^04^16^6, 5(2000) = 25°10'4r.'8. Tfie left panel shows observed 
(heavy histogram) and model (light smooth curve) spectra on a grid with 10" spacing along 
each axis. The velocity and temperature axis for each spectrum is illustrated in the lower 
left corner and runs from 6.3 to 8.0 km s~^in velocity and —0.5 to 3.3 K in temperature. The 
center panel shows contours of the total model N2H"'" column density (front and rear layers, all 
hyperfine components) in units of 10^'^ cm^^, overlaid on a grayscale of the integrated map of 
N2H+(101-012) emission from the combined map (lowest level and increment 0.2 Kkms~^). 
The right panel plots the velocity difference between the front and rear layers, Vi^, in units 
of m s"^. The crosses indicate those points for which there is sufficient optical depth and 
signal-to-noise to determine fin reliably. 
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TABLE 1 



Model Parameters 



Parameter 


Value 




5^04"^16'6 




25°11'06'.'8 


e 


25° 




12.0 K 


(TnT 


0.085 km 


X 


1 X 10-9 


?^pk,f 


1 X 10^ cm-3 




42000 AU 


^pk,r 


4 X 10^ cm"-^ 


/S.Zr 


4800 AU 


Ax 


2700 AU 


Ay 


1500 AU 


vo 


7.235 km s'^ 


dv/dx 


3.8 km pc~^ 



